library(readr)
file=read_delim("/mnt/md1200/6/zhaocunyou/Genome_analysis_by_wzj/All_sample_merged.HC.VQSR.genotype.heter.txt",col_names=T,delim="\t")
 ID1=c("X1T","X2B","SD1","SD1T","F15","F15T2","FC1","FC1T","FC6","FC6T","ZS6","ZS6T","JYS2","JYS2T","JYS1","JYS1T","FC2","FC2T","F14","F14T","FC3","FC3T","F10","F10T","F16","F16T","FC4","FC4T","CT1A","CT1B","CT2A","CT2B","CT3A","CT3B","CT5A","CT5B","JYS3","JYS3T","JYS4","JYS4T")
 ID2=c("X1T","X2B","M7","M8","M5","M6","M1","M2","M41","M42","M43","M44","M47","M48","M3","M4","M11","M12","M27","M28","M29","M30","M25","M26","M36","M35","M32","M31","M17","M18","M19","M20","M21","M22","M39","M40","M49","M50","M51","M52")
 
 disease=c(rep("SZ",4),rep("BD",10),"AD","AD","MDD","MDD","SZ","SZ",rep("BD",6),"AD","AD",rep("NC",8),rep("SZ",4))
 status=c(rep(c("DCD","DCC"),9),rep("CC",10),rep("NC",8),rep(c("DCD","DCC"),2))
 
 sample_info=data.frame(ID1,ID2,disease,status)
group1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39")
group1=unlist(strsplit(group1,"_"))
sample.info=sample_info[sample_info$ID2 %in% group1,]
 sample.id=as.character(sample.info$ID1)
 file1=file[,c("ID",names(file)[names(file) %in% sample.id])]

check=function(x){
 str1=unlist(strsplit(as.character(x[2])," "))
 if(str1[1]==str1[2]){ 
     return ("FALSE")
     }else{return ("TRUE")}
 }
 
 unitID=c()
 for(j in 2:15){
 test=data.frame(ID=file1$ID,file1[,j])
 test$group.heter=apply(test,1,check)
 test=tidyr::separate(test,names(test)[2],into=c("Ref","Alt"),sep=" ")
 test$unitID=paste(test$ID,test$Ref,test$Alt,sep=":")
 unitID=unique(c(unitID,as.character(test[test$group.heter=="TRUE",]$unitID)))
 }
 length(unitID)
 
 id=data.frame(id=unitID)
write.table(id,"/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/all.heter.genotype.snp",quote=F,row.names=F,sep="\t")
id=read.table("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/all.heter.genotype.snp",head=T,sep="\t")

id2=tidyr::separate(id,id,into=c("Chr","Pos","Ref","Alt"),sep=":")
id2$id=paste(id2$Chr,id2$Pos,sep=":")
write.table(id2,"/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/all.heter.genotype.snp",quote=F,row.names=F,sep="\t")

con=id2
filea=read.csv("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/20201207/all.FDR.sig.at.least.one.add.direction.same.diff.csv",head=T)
filea$id=paste(filea$Chr,filea$Start,sep = ":")
con=con[!con$id %in% as.character(filea$id),]
conid=as.character(id2$id)
conid.random=sample(conid,150000,replace = FALSE)
con1=con[con$id %in% conid.random,]
con1=data.frame(con1$Chr,con1$Pos,con1$Pos,con1$Ref,con1$Alt)
write.table(con1,"/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/20201207/con_genotype.random1.txt",quote=F,row.names=F,sep="\t")
